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We present the results of the physical point simulation in 2+1 flavor lattice QCD with the non- 
perturbatively 0(a)-improved Wilson quark action and the Iwasaki gauge action at ^ = 1.9 on a 
32 3 x 64 lattice. The physical quark masses together with the lattice spacing is determined with 
m w , tuk and ma as physical inputs. There are two key algorithmic ingredients to make possible the 
direct simulation at the physical point: One is the mass-preconditioned domain-decomposed HMC 
algorithm to reduce the computational cost. The other is the reweighting technique to adjust the 
hopping parameters exactly to the physical point. The physics results include the hadron spectrum, 
the quark masses and the pseudoscalar meson decay constants. The renormalization factors are 
nonperturbatively evaluated with the Schrodinger functional method. The results are compared 
with the previous ones obtained by the chiral extrapolation method. 
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I. INTRODUCTION 



The physical point simulation is one of the essential 
ingredients in the first principle calculation of lattice 
QCD. However, it is still a tough challenge because of 
the rapid growth of the computational cost with the up- 
down (ud) quark mass reduced toward its physical value. 
At present simulation points are typically restricted to 
m,r>250 MeV. The most popular strategy to obtain the 
results at the physical point is chiral extrapolation with 
the use of chiral perturbation theory (ChPT) as a guiding 
principle. This strategy, however, has several problems: 
(i) It is numerically difficult to precisely trace the loga- 
rithmic quark mass dependence of the physical quantities 
predicted by ChPT. (ii) It may not be always possible to 
resort to ChPT as a good guiding principle for chiral ex- 
trapolation, (iii) The kinematics changes as the quark 
mass increases. A typical example is the p — > tttt de- 
cay which is not allowed for the increased ud quark mass 
away from the physical value, (iv) Our final destina- 
tion is to incorporate the different up and down quark 
masses. The isospin breaking effects are so tiny that re- 
liable evaluation would be difficult by the chiral extrap- 
olation method. 

In this article we present the results of the phys- 
ical point simulation which has been pursued as the 
PACS-CS project based on the PACS-CS (parallel ar- 
ray computer system for computational sciences) com- 
puter with a peak speed of 14.3 Tflops developed at 



University of Tsukuba [l|-Q. The simulation is carried 
out with the nonperturbatively 0(a)-improved Wilson 
quark action^ and the Iwasaki gauge action^ on a 
(3 fm) 3 box at the lattice spacing of a = 0.08995(40) 
fm. There are two types of problems in the physical 
point simulation. First, we need to reduce the compu- 
tational cost which rapidly increases as the ud quark 
mass decreases. This difficulty is overcome thanks to 
the domain-decomposed HMC (DDHMC) algorithm |q| 
with the mass-preconditioning[3, @. In Refs. 0, [l(| 
this algorithm was successfully applied to investigate 
the chiral behaviors of the pseudoscalar meson sector 
and the hadron masses including both the mesons and 
the baryons, where the pion mass covers from 156 MeV 
to 702 MeV. The second problem is fine-tuning of the 
quark masses to the physical point after we reach around 
the physical point. This task is accomplished with the 
reweighting technique which allows us to cover a small 
variation of simulation parameters in a single Monte 
Carlo run[ll|. We explain the details of the method and 
present the physics results on the physical point without 
interpolation or extrapolation. 



*Present address: Theoretical Physics Laboratory, RIKEN, Wako 
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This paper is organized as follows. In Sec. |TT] we 
present the simulation details including the parameters 
and the algorithm. Section IIIII is devoted to describe the 
reweighting method. We present the physics results on 
the physical point in Sec. IIVI Our conclusions are sum- 
marized in Sec. W\ 
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II. SIMULATION DETAILS 
A. Actions 

We employ the Iwasaki gauge action Q and the non- 
perturbatively Q(a)-improvcd Wilson quark action as in 
the previous works [3, 13 ■ The former is composed of a 
plaquette and a 1 x 2 rectangle loop: 



with a = -0.331 and c = 1 - 8c x = 3.648. The latter 
is expressed as 



5 g = i| c ° E trt/ > 

I plaquette 



P l + Ci 2^ trU rtg } (!) 

rectangle 
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(2) 



where we consider the case of a degenerate up and down 
quark mass k u = k^. The Euclidean gamma matrices 
are defined in terms of the Minkowski matrices in the 
Bjorken-Drell convention: jj = —i"f J BD (j = 1, 2, 3), 74 = 

1b Di 75 = Ibd and °> = WlvtilA- Tne neld strength 
F^u in the clover term is given by 

, 4 

= iE^(^( n )-^K)> ( 3 ) 

1=1 

C/i(n) = U^Un+^uU^^, (4) 
C/ 2 (n) = l^X-^X-A,^"-^ (5) 

= ul_ Pi<ll ul_ ji _ ^u n -p,-{ > ^u n - i)!V , (6) 

C/ 4 (n) = ul^U^Un+^U^. (7) 

The improvement coefficient csw for 0(a) improvement 
was determined nonperturbatively in Ref. Q. 

B. Simulation parameters 

Simulations are performed employing the same pa- 
rameters as in the previous workQ: a 32 3 x 64 lat- 
tice at = 1.90 with c sw = 1-715 Q. We choose 
{k u<1 ,k s ) = (0.137 785 00,0.136 600 00) for a degenerate 
pair of up and down quarks and a strange quark. This 
combination of the hopping parameters was supposed to 
be the physical point based on the analysis of the pre- 
vious results @. The lattice spacing is determined as 
a = 0.08995(40) fm from the m^, tuk , mn results on the 
physical point after the reweighting procedure. Table U 
summarizes the run parameters. After thermalization we 
calculate hadronic observables solving quark propagators 



at every 20 trajectories (5 MD time units), while we mea- 
sure the plaquette expectation value at every trajectory. 
The reweighting factors for the up-down and the strange 
quarks are evaluated at every 100 trajectories (25 MD 
time units). The choice of sparse measurements is due 
to the demanding computational cost of the reweighting 
factors. The hadronic observables measured at every 20 
and 100 trajectories show consistency within error bars. 
For the pion mass we find that the former has larger mag- 
nitude of error: 0.0719(37) and 0.0693(27). This could 
be due to wavy behavior of pion propagators on a couple 
of configurations caused by the statistical fluctuation. 



C. Algorithm 

Our base algorithm for the degenerate up-down quarks 
is the DDHMC algorithm^ which makes a geometric 
separation of the up-down quark determinant into the 
UV and the IR parts with the domain-decomposition of 
the full lattice into small blocks. This UV/IR separa- 
tion naturally introduces the multiple time integration 
scheme[l3j in the molecular dynamics (MD) steps. We 
employ the nested simple leapfrog with QPQ ordering for 
the multiple time step MD integrator. According to the 
relative magnitude of the force terms coming from the 
gauge part and the UV and the IR parts of the up-down 
quarks we choose the associated step sizes such that 

5r g ||F g || « ^uvll-Fuvll « 5r IR ||F IR ||, (8) 

where 5r g = t/NqN^, St vv = t/N^, St ir = t/N 2 
with r the trajectory length and (No, N\, N2) a set of 
integers to control the step sizes. 
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In the previous work we used the mass-preconditioned 
DDHMC (MPDDHMC) algorithm for the run at 
(n ud ,n s ) = (0.137 810 00,0.136 400 00) which gives the 
lightest up-down quark mass@. A preconditioner con- 
trolled by an additional hopping parameter K' ud = pK u d 
is incorporated to tame the fluctuation of the IR force Fir 
in the original DDHMC algorithm by dividing it into Fir 
and F IR . The former is derived from the preconditioned 
action and the latter from the preconditioner. In this 
work we employ twofold-mass-preconditioned DDHMC 



ll^gll : ll^uvll : KrII : KrII 



with p! = 0.9995 and p 2 = 0.9900. We 
choose (N ,Ni,N 2 ,N 3 ,N^) = (4,4,2,4,4) for the 
associated step sizes: Sr g = T/N0N1N2N3N4, 
6t vv = r/N^NaNi, 5t& = t/N 2 N z N a , St{ r = 
t/A^A^jJtir = t/N±. This choice results in rather high 
acceptance rate found in Table Q] The replay trick [(I [l4| 
is not incorporated. 

For the inversion of the Wilson-Dirac operator dur- 
ing the MD steps we implement the same algo- 
rithmic techniques as for the run at (k uc [,k s ) = 
(0.13781000,0.13640000) in the previous work@. 
There are three important points to be noted. First, 
the initial solution vector is provided by the chronolog- 
ical guess with the last 16 solutions [Hj]. We demand a 
stringent stopping condition \Dx — b\/\b\ < 10~ 14 to as- 
sure the reversibility. Second, the inversion algorithm 
is a nested BiCGStab solver consisting of an inner and 
an outer solvers. The former plays the role of a precon- 
ditioner whose calculation is accelerated by single pre- 
cision arithmetic with an automatic tolerance control 
ranging from 10~ 3 to 10~ 6 . The latter is implemented 
with double precision imposing a stringent tolerance of 
10~ 14 . Third, the deflation technique is incorporated in 
a nested BiCGStab algorithm: Once the inner solver be- 
comes stagnant during the inversion of the Wilson-Dirac 
operator, the solver algorithm is automatically replaced 
by the GCRO-DR (generalized conjugate residual with 
implicit inner orthogonalization and deflated restarting) 
algorithm flU]. This saves us from the difficulties due 



(MP 2 DDHMC) algorithm which split F m into F m , F/ R 
and F IR . This decomposition is controlled by two addi- 
tional hopping parameters n[ ld = p\K and /t" d = p\p 2 n. 
Fir is derived from the action preconditioned with K' ud . 
The ratio of two preconditioned with n' ud and ft" d gives 
Fj' R . Fjjj is from the heaviest preconditioners with K^ d . 
We find the following relative magnitude for the force 
terms: 



:||Fir||« 16:4: 1:1/7: 1/60 (9) 



I 

to possible small eigenvalues allowed in the Wilson-type 
quark action. We refer to Appendix B in Ref. @ for more 
details of the inversion algorithm. 

The strange quark is simulated with the UV-filtered 
PHMC (UVPHMC) algorithm[l3-[20i where the action 
is UV-filtered[ll| after the even-odd site preconditioning 
without domain-decomposition. We set the step size as 
St s = <5r IR according to our observation ||F S || w ||Fir||. 
This algorithm is made exact by correcting the polyno- 
mial approximation with the global Metropolis testj22j 
at the end of each trajectory. In Table Q] we find that the 
choice of Af po i y = 220 yields 95% acceptance rate. 



III. RE WEIGHTING METHOD 
A. Formalism 

Let us consider evaluating (0[[/](«* d , K s))« d ,Kj), 
which is the expectation value of a physical observable 
O at the target hopping parameters («:* d ,K*), using the 
configuration samples generated at the original hopping 
parameters (K u d, n s ). We assume that p u d = ft u d/K* d — 1 
and p s = k s /k* ~ 1. With this assumption, the expec- 
tation value is rewritten as follows using the single his- 
togram rcweighting method[ll|: 
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/Pt/0[[/]« d ,<)|det[^ <d [[/]]| 2 det[^[[/]]e- g 8 ^ 
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) (K ud ,K 


D K * [U] 
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|det[^ ud [C/]]pdct[^ Ks [C/]]e- s 9 ^ 
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( J R ud [C/]i? s [C/]) (Kud , Ks) 



where the reweighting factors arc defined as 



Rud[U] = 
R S [U] = 



det 



det 



d«'Ju] 

D K [U 



[U] 



and 



D Kq [U] = I + Kq {T + M) (g = ud,s) 



(11) 
(12) 

(13) 



with T the local clover term including the nonperturba- 
tive csw and M the hopping matrix. The above expres- 
sion (fTU)) demands us to evaluate the reweighting factors 
R u d[U] and R S [U] on each configuration. For later con- 
venience wc define 



W[U]( Pq ) EE 



D Kq [U] 



(14) 



With /9q = Kq/fi*. 



B. Evaluation of reweighting factors 



The ratio W 1 is further simplified as follows: 

D Kud P] 



= Pnd + (l- Pud )D-}[U] (17) 



with the use of D Kud [U] = pudD <d [U] + (1 - p ud ). We 
just need D~* to calculate W~ x . 

For the strange quark we assume that det [W[J7](p s )] 
is positive. The corresponding reweighting factor is eval- 
uated as 

R S [U] = det [W[U](p a )] 



JP7 ? tX> ?ye -|W" 1/2 [y](Ps)r,| 2 + |»)| 2 -|»)| : 



\W"^ 2 [u]( Pa ) v \ 2 +\n\ 2 



IV 



(18) 



The reweighting factor i? u d[C^] can be evaluated with a 
stochastic method. Introducing a complex bosonic field 
77, whose spin and color indices are suppressed here, the 
determinant of W is expressed as 

R nd [U] = |det[Fy[C/]( Pud )]| 2 



[VrfVr)e-\ w ~ 1 W(P" d '> r >\ 2+ \ r >\ 2 -\ r >\'" 
J VrfVr)e-W 2 



(15) 



means the expectation value with respect to 
, Nn) which arc random 



where (• 

T). Given a set of 77W (i = 1, 
noises generated according to the Gaussian distribution, 
the reweighting factor is evaluated as 



Rud[U] 



lim 



N„ 



1 N n 



|lV- 1 [C/](Pud)»?| 2 + l'7| : 



(16) 



With the assumption of p s ~ 1 we expect that W[C/](/9 s ) 
is so close to the identity matrix that its eigenvalues are 
enclosed by a unit circle centered at (1,0) in the complex 
plane. In this case we can evaluate W~ 1 ^ 2 [U](p s )r] by the 
Taylor expansion around identity. 

To evaluate the matrix square root W^ 1 ^ 2 [U}(p s ) we 
first parametrize W~ 1 [U](p s ) as 



w-'PKps) = 



DM 

D K ,[U] 

p s + (l- Ps )D-}[U] 
l-(l-p s ) (l-D-}[U] 
l-X[U}(p s ) 



(19) 



where |1 - p s \ -C 1 and ||X[77](p s )|| < 1. We em- 
ploy the recursive expression for the Taylor expansion 

of w-v^iriips)^: 
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CO 



c ±x 

CO 



Cl 



C2 



r\ -\ X 

CN-2 



Cjy 
cn-1 



Xr] 



(20) 



where the argument [£/](p s ) for the matrices is sup- 
pressed. The coefficients are given by 



2.7 



(21) 



with Co = 1. The advantage of the recursive procedure 
is to reduce the round-off errors in the summation from 
the lower-order to the higher-order contributions in the 
Taylor expansion. The truncation error and the order 
of the Taylor expansion N arc monitored and controlled 



during the simulation by explicitly evaluating the resid- 
ual r = WiW^^W- 1 / 2 - W- 1 )!]]]/ /\\r]\\. We enforce the 
condition r < 10~ 14 for N. 



To reduce the fluctuations in the stochastic evalua- 
tion of i? u d[C/] and R S \U] we employ the determinant 
breakup technique [HI, [5"" 
K* is divided into Nb subintervals 



The interval between n q and 



K g + (N B - 1)A„ k*} with A q = (k* - Kq )/N B . Thus 
the determinant of iy[[/](p g ) is broken up as 



det [W[U](p q )] = det 



W[U] 



det 



W[U] 



2A C 



■det 



W[U] 



+ (N B - 1)A 9 



(22) 



where each determinant in the right hand side is eval- 
uated with an independent noise set of r\. For strange 
quark reweighting, W~ x are simply replaced by W -1 ! 2 
in Eq. ([22]). 



C. Parameters and results for reweighting factors 

Our choice of the target hopping parameters are 
«d><) = (0.137 796 25,0.136 633 75). The subintervals 
for the determinant breakup are A u d = (0.137 796 25 — 
0.137 785 00)/iV B with N B = 3 for the ud quark and 
A s = (0.136633 75- 0.136 600 00)/N B with N B = 3 for 
the s quark. Each piece of the divided determinant is 
evaluated stochastically employing 10 sets of r/ at every 
100 trajectories (25 MD time units). The order of Taylor 
expansion N was mostly 5 for each of the strange quark 
reweighting break up. 

Figure [1] shows configuration dependence 
of the reweighting factors from (ft u d, k s ) = 
(0.137 785 00,0.136 600 00) to « d X) 
(0.137 796 25,0.136 633 75) which are normalized as 
(Rud,s) = 1 an d (R u dRs) = 1- The fluctuations of R u d 
and Rg are within a factor of 10. Their product has 
slightly amplified fluctuations. In Fig. [2] we plot the 
reweighting factors as a function of the plaqucttc value on 
each configuration. An important observation is a clear 



correlation between the reweighting factors and the pla- 
quette value: The former increases as the latter becomes 
larger. Thanks to this correlation the distribution of the 
plaquette value at « d ,K*) = (0.137 796 25,0.136 633 75) 
is moved in the positive direction. This is the expected 
behavior, because the target hopping parameters arc 
larger than the original ones. The situation is quan- 
titatively illustrated in Fig. [3l where the reweighted 
plaquette values with i? u( j and R s arc individually 
plotted as a function of the number of noise. The results 
look converged once the number of noise goes beyond 
four. 

Since the formula of Eq.[TU]is the identity, the reweight- 
ing procedure is always assured if we have infinite statis- 
tics. In case of finite statistics in practical simulations, 
however, we should be concerned with the possible situ- 
ation that the original and the target points are far away 
such that the distributions of observable fail to overlap 
each other. This problematic case could be detected by 
monitoring the behavior of the expectation value for the 
observable as the reweighting parameters are monoton- 
ically moved from the original point: The expectation 
value of the observable stops varying with diminishing 
error bar. To check the reliability of our reweighting 
procedure we have investigated the behavior of the ex- 
pectation value of the plaquette against the reweighting 
with respect to the strange quark from k s = 0.136 600 00 
to k s = 0.136 690 00 with N B = 4 and 8, the latter of 
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which yields the same amount of breakup interval as 
A s = (0.136 633 75 - 0.136 600 00)/3 = 0.000 01125 in 
our choice. Since the plaquette value has much narrower 
distribution than the hadron propagators at each time 
slice, this is a stringent test to check the overlap of the 
distributions of the observable at the original and the tar- 
get points. Figured] shows the behavior of the reweightcd 
plaquette value evaluated with 10 noise sources as a func- 
tion of the reweighting parameter k s . We do not observe 
any sign that the reweighted plaquette value stagnates 
at some point: It shows almost linear behavior with con- 
stant magnitude of error up to k s = 0.136 690 00 which 
is far beyond the physical point of k s = 0.136 633 75. 
Furthermore Nb = 4 and 8 cases give consistent re- 
sults. In Fig. [5] we plot the reweighting factor R s from 
k s = 0.136 600 00 to k s = 0.136 690 00 with N B = 4 and 
8 as a function of the plaquette value on each configura- 
tion, which is normalized as (R s ) = 1. Both cases show 
quite similar distributions, which confirm that our choice 
of breakup interval A s = 0.000 011 25 is sufficiently small. 
In Fig. [6] we also present the reweighted plaquette value 
with R s as a function of the number of noise. The results 
with Nb = 4 and 8 become fairly consistent once we em- 
ploy more than two noise sources. We have repeated the 
same analyses for the reweighting with respect to the ud- 
down quark from k s = 0.137 785 00 to k s = 0.137 800 00 
with Nb = 2 and 4. The same conclusion is obtained as 
in the strange quark case. This is easily expected from 
similar behaviors for i? u d and R s found in Figs. [1] [2] and 

El 



IV. HADRONIC OBSERVABLES 

A. Hadron masses, quark masses and decay 
constants at simulation point 

We measure the meson and the baryon correlators em- 
ploying appropriate operators. The general form of the 
meson operators is expressed as 



Ml 9 (x)=q f (x)Tq g (x), 



(23) 



where / and g denote quark flavors and T are 16 Dirac 
matrices T = I, 75, 7^, 17^75 and i[y fl ,^ l ,]/2 {p,, v = 
1, 2, 3, 4). The octet baryon operators are given by 

Of a « h (x) = e abc {{q a f {x)) T C lb q b {x))ql a {x), (24) 

where a, b, c are color indices, C — 7472 is the charge 
conjugation matrix and a, = 1, 2 labels the z-componcnt 
of the spin 1/2. The S- and A- like octet baryons are 
distinguished by the flavor structures: 



E-like 
A-like 



Q[fh]g _|_ Q[gh]f 
Q[fh\g _ Q[gh]f _ 20lfg] h 

' 



(25) 



(26) 



where 0^ h = O fgh - O afh . We define the decuplet 
baryon operators for the four z-components of the spin 
3/2 as 



ttqWfcr+ftxM^x), (27) 

hi 



D{%(x) = e abc [{{q a f {x)) T CT Q q b g {x))ql x {x) 

-{{q%x)) T CT + q b g {x))ql 2 {x)]/Z, (28) 
Dtf /2 (x) = e abc [{{q a f (x)) T Cr q b g {x))q c h2 {x) 

-((q a f (x)) T Cr_q b g (x))q c hl (x)]/3, (29) 
e abc ((q](x)) T CT. q b (x))ql 2 (x), (30) 



D%{x) 



where r± = (7i=F72)/2, To = 73 and the flavor structures 
should be symmetrized. 

The meson and the baryon correlators are calculated 
with point and smeared sources and a local sink. The 
smeared source is constructed with an exponential smear- 
ing function vE , (|a?|) = A q exp(—B q \x\) (q = ud,s) where 
\l/(0) = 1 for the ud and s quark propagators. Employ- 
ing a couple of thermalizcd configurations we adjust the 
parameters such that the pseudoscalar meson effective 
masses reach a plateau as soon as possible. Our choice is 
^ ud = 1.2, B ud = 0.07 and A s = 1.2, B s = 0.18. 

To reduce the statistical error of the zero mo- 
mentum hadron correlators we employ two meth- 
ods. One is the choice of four source points at 
(x ,y ,zo,io)=(17,17,17,l), (1,1,1,9), (25,25,25,17), 
and (9,9,9,25). The other is the use of possible spin 
states: three polarization states for the vector meson and 
two (four) spin states for the octet (decuplet) baryons. 
The correlators with different sources and spin states are 
averaged on each configuration before the jackknife anal- 
ysis. 

Figure [7] shows effective mass plots for the meson and 
baryon propagators with the smeared source, where we 
assume a single hyperbolic cosine function for the for- 
mer and a single exponential form for the latter. We ob- 
serve good plateaux starting at small values of t, showing 
that the excited state contributions are suppressed. The 
hadron masses are extracted by uncorrelated \ 2 fit to the 
propagators, since we find instabilities in correlated fit us- 
ing covariance matrix. The horizontal bars in Fig. [7] rep- 
resent the fit ranges, which are [t m in,^max] = [13,30] for 
the pseudoscalar mesons, [10, 20] for the vector mesons 
and [6, 10] for the baryons, and the resulting hadron 
masses with 1 standard deviation error band. The nu- 
merical values are summarized in Table [XT] The statis- 
tical errors are estimated with the jackknife method. In 
Fig. [5] we plot the binsize dependence of the error for 
the pseudoscalar meson masses. The magnitude of error 
shows flat behaviors against the binsize within the er- 
ror bars. Since similar binsize dependences are found for 
other particle types, we employ a binsize of 100 MD time 
(4 gauge configurations) for the jackknife analysis. As a 
cross check we also carry out the bootstrap error estima- 
tion with 5000 samples. For all the physical quantities at 
the original and the target points the bootstrap samples 
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show clear normal distribution and the error estimates 
agree with those of the jackknife method within 10%. 

For the quark masses and the decay constants we have 
accomplished an important improvement since the pre- 
vious publication]^]: a nonperturbative determination of 
renormalization factors based on the Schodinger func- 
tional scheme [IB-US]- The bare quantities are calculated 
with the same method as in Ref. @. 

The bare quark mass is defined by the axial vector 
Ward-Takahashi identity (AWI): 



_ AWI , - AWI 

rrif + nig 



|V 4 Af p |PS) 
(0\P\PS) 



(31) 



where P is the pseudoscalar operator and |PS) denotes 
the pseudoscalar meson state at rest consisting of / 
and g (/, g = ud, s) valence quarks. The axial vector 
current is nonperturbatively 0(a)-improved as A™ p = 
A± + C4V4P with V4 the symmetric lattice derivative 
and c A = -0.038 761 060. The ratio of the matrix ele- 
ments is evaluated by 



fhj wl + m£ WI 



m PS 



Ci 



(32) 



where mps, C S A and Cp are extracted from a simultane- 
ous x 2 n t to 

S inh(-m PS (t-r/2)) 



exp(mpsT/2) 



and 



(P(t)P s {0)) = 2C] 



cosh(-TO PS (t - T/2)) 
exp(TO PS T/2) 



(34) 



with P s the smeared pseudoscalar operator and T = 64 
the temporal extent of the lattice. The fit ranges are 
chosen to be [t m i n , i max ] = [13,25] for the former and 
[13, 30] for the latter. The renormalized quark mass in 
the continuum MS scheme is defined as 



MS _ 7MS-AWI 
m f ~ A m m f > 



(35) 



where Z^ s = Za/Zp is nonperturbatively determined 
in the Schrodinger functional scheme. In Table [IV 
we present the results for mjjjf and to^ s renormal- 
ized at fi = 2 GeV together with the correspond- 
ing bare quark masses to^ wi and mf WI . We use 

ZjP(2GeV)=1.441(15)|22|. The statistical errors are es- 
timated by the jackknife analysis with the choice of the 
same binsize as for the hadron masses. 

The bare pseudoscalar meson decay constant defined 
by 



(0|A™ p |PS) 



rbarc 
J PS 



m PS . 



is evaluated from the following formula: 



rbarc 
J PS 



C\ 



CI 



2 Cl 



TOPS 



(36) 



(37) 



We extract tops, C% , Cp and C l P from a simultaneous 
fit of Eqs. ([33L ([341 and 



(P(t)P l (0)) = 2C 



l cosh(-m PS (i-T/2)) 
exp(m PS r/2) 



(38) 



with P l the local pseudoscalar operator. The fit ranges 
are [13,25], [13,30] and [15,25], respectively. The renor- 
malization is given by 



r ry /-bare 

IPS — ^A/ps > 



(39) 



with Za = 0.8563(52) [23 the nonperturbative renormal- 
ization factor in the Schrodinger functional scheme. In 
Table [TV] we list the results for / PS and /^| rc with the 
statistical errors evaluated in the same manner as for the 
quark masses. 



B. Hadron masses, quark masses and decay 
constants at target point 

In Fig. [5] we present the effective masses for the 
reweighted meson and baryon propagators with the 
smeared source. Comparing with the original case in 
Fig.[7]the error bars are slightly enlarged by the reweight- 
ing procedure. We apply the uncorrelated \ 2 fit to the 
reweighted hadron propagators at the target point choos- 
ing the same fit ranges and jackknife binsize as in the sim- 
ulation point. The results are summarized in Tables [IT] 
and IIIIl where we also present the previous results ob- 
tained by the chiral extrapolation method in Ref. Q for 
comparison. 

To investigate the reweighting effects on the hadron 
effective masses, we show the effective masses for the 
pseudoscalar mesons with and without the reweighting 
factors in Fig. 1 1 01 where ?7 SS is a fictitious pseudoscalar 
meson consisting of two strange quarks. For all the cases 
the partially quenched results (PQ) show lighter effective 
masses than the unitary results at the simulation point. 
They are further reduced by the reweighting procedure 
(PQ+RW). For other hadron channels the reweighting 
effects are less clear partly because of the larger error 
bars. 

In Fig. [TT] we plot the n, p and nucleon masses at the 
target point as a function of the number of noise. The 
situation is quite similar to the plaquettc case: Five or 
six noises appear sufficient to obtain a reliable estimate. 
This is also the case for other hadron masses. 

Figure [T^] compares the measured hadron masses nor- 
malized by mo with the experimental values. The re- 
sults for TOtt/toq and rriK/mn, which are sizably devi- 
ated from the experimental values at the simulation point 
(black symbols), are properly tuned to the physical val- 
ues within error bars at the target point of (K* d ,«*) = 
(0.137 796 25,0.136 633 75). The lattice spacing is deter- 
mined as a — 0.08995(40) fm from toq. A large dis- 
crepancy found for m p /mn may be resolved by a proper 
treatment of p meson as the resonance [10, Hl|. We plan 
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to do so for the p, K* mesons and A baryon. For other 
hadron masses we find less than 5% deviation from the 
experimental values. An increasingly larger deviation ob- 
served for lighter baryons may be due to finite size effects. 

Possible finite size effects on the pseudoscalar meson 
masses based on the NLO formulae of ChPT[32[ are dis- 
cussed in Sec. IV D of Ref. Q. The expected correc- 
tions are less than 2% for and m K at the physical 
point. The magnitude is smaller than the statistical er- 
rors found in Table IIIIl For the baryon masses the heavy 
baryon ChPT predicts less than 1% corrections at the 
physical point on our physical volume as listed in Table 
X of Ref. 0. 

Although Fig. [12] clearly shows that further tun- 
ing is not really necessary, it would be instructive 
to pin down the physical point in the (1/k uc i, 1/k s ) 
plane. The physical point plotted in Fig. [13] is de- 
termined by a combined linear fit of (m^/ma) 2 and 
(m K /m n ) 2 at (K ud ,K s ) = (0.137 785 00,0.136 600 00), 
(K* d ,/c*) = (0.137 796 25,0.136 633 75) and two more 
reweighted points given by (0.137 796 25, 0.136 633 75 ± 
A s ). The fit functions are 



/ \ 2 K K 

frnx\ = c K + ± + 2^ 
\ ma, J n ud k s 



(40) 



(41) 



with 2 free parameters. The experimental values 
of m^/mn and mx/mri are reproduced at (K n d,K s ) — 
(0.137797(4), 0.136 635(16)), whose central value is al- 
most exactly hit by our target point (k*j,k*) = 
(0.137 796 25,0.136 633 75). 

The quark masses and the pseudoscalar decay con- 
stants are extracted by repeating the same analyses 
as in the simulation point. The results are sum- 
marized in Table IIVI The quark masses are deter- 
mined as m^f (2 GcV)=2.97(28)(03) MeV and mp(2 
GeV)=92.75(58)(95) MeV with a" 1 = 2.194(10) GeV, 
where the second error is due to the nonperturbative 
renormahzation factor obtained by the Schrodinger func- 
tional method (27| . We find that our quark masses are 
comparable to recent estimates in the literature [33}. The 
discrepancy between the quark masses in this work and 
those in Ref. [9] is mainly due to the difference in 
the renormahzation factors. The nonperturbative esti- 
mate gives about 30% larger value than the perturba- 
tive one[27J. For the pseudoscalar meson decay con- 
stants we obtain / w = 124.1(8.5)(0.8) MeV and f K = 
165.5(3.4)(1.0) MeV with the second error coming from 
the nonperturbative renormahzation factor [27| . These 
values should be compared with experiment: f v = 
130.4 ± 0.04 ± 0.2 MeV and f K = 155.5 ± 0.2 ± 0.8 ± 0.2 
MeV[29|. Note that the NL0 ChPT analyses predict 4% 
(1.5%) deficit for f v (fx) on a (3 fm) 3 box at the physical 
point due to the finite size effects [9l. |32|. 

V. CONCLUSION 

We have presented the results of the physical point 
simulation in 2+1 flavor lattice QCD with the 0{a)- 



improvcd Wilson quark action. This is accomplished 
by two algorithmic ingredients: the DDHMC algorithm 
with several improvements and the reweighting tech- 
nique. The former contributes to cost reduction and the 
latter is required for fine-tuning to the physical point. 

Clear reweighting effects arc observed on several ob- 
serbablcs: The plaquette value increases and the hadron 
masses are reduced in agreement with the expecta- 
tion for the reweighting from the simulation point at 
(«ud, ^s) = (0.137 785 00, 0.136 600 00) to the target point 
at KdX) = (0.137 796 25,0.136 633 75). We are al- 
lowed to properly tune the measured values of m n , rriK 
and ma to their experimental ones. 

We extract the hadron masses, the quark masses and 
the pseudoscalar decay constants directly on the physi- 
cal point after the reweighting procedure. For the hadron 
masses we find less than 5% deviation from the experi- 
mental values except the p meson case which requires 
a proper analysis as the resonance. The results for the 
quark masses renormalized in the MS scheme at p = 2 
GeV are presented with the nonperturbative renormal- 
ization factor determined by the Schrodinger functional 
method. The large enhancement of the quark masses 
compared to those in Ref. @ is attributed to the differ- 
ence between the nonperturbative renormahzation factor 
and the perturbative one. 

The physical point simulation, which has been the 
long-standing problem in lattice QCD, is achieved in this 
work. It appears to us that it is not worthwhile to in- 
crease the statistics with the present simulation param- 
eters. More important as the next step is to repeat the 
physical point simulation with larger and finer lattices. 
Further reduction of the finite size effects and the finite 
cutoff effects will make possible precision measurements 
of physical observables at 1% level. 
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TABLE I: Simulation parameters. MD time is the number of trajectories multiplied by the trajectory length r. 





0.137 785 




0.136 600 


#run 


5 


T 


0.25 


block size 


8 4 


(N ,Ni,N2,N 3 ,m) 


(4,4,2,4,4) 


Pi 


0.9995 


Pi 


0.9900 


iVpoly 


220 


Replay 


off 


MD time 


2000 


(P) 


0.571 082(9) 


(e- dH ) 


0.9916(81) 


Pacc(HMC) 


0.8109(45) 


Pacc(GMP) 


0.9519(27) 



TABLE II: Meson and baryon masses in lattice units at original and target points. 





original 


target 


physical point in Ref. [9J 




0.137 785 00 


0.137 796 25 






0.136 600 00 


0.136 633 75 




7T 


0.0693(27) 


0.0617(28) 


0.0620(9) 


K 


0.2321(10) 


0.2270(9) 


0.2287(33) 




0.3203(7) 


0.3138(6) 


0.3168(43) 


p 


0.331(38) 


0.272(39) 


0.357(16) 


K* 


0.4028(55) 


0.393(11) 


0.4118(72) 


4> 


0.4652(17) 


0.4605(28) 


0.4634(61) 


N 


0.441(12) 


0.447(21) 


0.438(20) 


A 


0.5147(63) 


0.518(10) 


0.502(10) 


E 


0.5485(38) 


0.5484(62) 


0.531(11) 




0.6022(27) 


0.6001(28) 


0.5991(75) 


A 


0.593(16) 


0.587(27) 


0.587(19) 


E* 


0.6557(67) 


0.658(12) 


0.657(15) 




0.7114(39) 


0.7113(53) 


0.718(12) 


n 


0.7655(34) 


0.7624(34) 


0.769(11) 
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TABLE III: Meson and baryon masses in physical units at target point. Experimental value for m?, BS is estimated by m< 
\/2m 2 K -m|. 





target [GeV] 


physical point in Ref. [9] [GeV] 


experiment [GeV] [29] 




0.137 796 25 








0.136 633 75 






7T 


0.1354(62) 




0.1350 


A 


0.4980(22) 




0.4976 




0.6884(32) 


0.6895(20) 


0.6906 


P 


0.597(86) 


0.776(34) 


0.7755 


A" 


0.861(23) 


0.896(9) 


0.8960 


<t> 


1.0102(77) 


1.0084(40) 


1.0195 


N 


0.982(45) 


0.953(41) 


0.9396 


A 


1.137(25) 


1.092(20) 


1.1157 


E 


1.203(11) 


1.156(17) 


1.1926 




1.3165(60) 


1.304(10) 


1.3148 


A 


1.289(59) 


1.275(39) 


1.232 


E* 


1.444(25) 


1.430(23) 


1.3837 




1.560(10) 


1.562(9) 


1.5318 


Q 






1.6725 



TABLE IV: Quark masses and pseudoscalar decay constants at original and target points. Renormalization factors are non- 
perturbative in this work, while perturbative in Ref. • 





original 


target 


physical point in Ref. [9] 


experiment [291 




0.137 785 00 


0.137 796 25 








0.136 600 00 


0.136 633 75 








0.001241(95) 


0.000 939(87) 


0.001042(32) 




arh* WI 


0.030 44(9) 


0.029 34(12) 


0.029 99(70) 




mfi [MeV] 


3.92(30) (04) 


2.97(28)(03) 


2.527(47) 




mf s [MeV] 


96.23(52)(98) 


92.75(58)(95) 


72.72(78) 




?n s /m ud 


24.5(1.8) 


31.2(2.7) 


28.78(40) 




/■bare 
"At 


0.0701(35) 


0.0661(45) 


0.0753(22) 






0.0898(16) 


0.0881(19) 


0.0897(18) 




U [MeV] 


131.7(6.6)(0.8) 


124.1(8.5)(0.8) 


134.0 (4.2) 


130.4 ±0.04 ±0.2 


fx [MeV] 


168.7(2.7) (1.0) 


165.5(3.4)(1.0) 


159.4(3.1) 


155. 5 ±0.2 ±0.8 ±0.2 


.fi</U 


1.280(60) 


1.333(72) 


1.189(20) 
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FIG. 1: Configuration dependence of reweighting factors from (K u d,K s ) = (0.137 785 00,0.136 600 00) to (s'd, K s) 
(0.137 796 25,0.136 633 75). 
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FIG. 2: Reweighting factors from (n ud ,n s ) = (0.137 785 00,0.136 600 00) to («£d. «S) = (0.137 796 25,0.136 633 75) as a function 
of plaquette value. 
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FIG. 3: Reweighted plaquette values with R u d and R B as a function of the number of noise. 
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FIG. 4: Reweighted plaquette values with R a as a function of target value of k s . Interval from k s = 0.136 600 00 to 0.136 690 00 
is divided by Nb = 4 (black) and 8 (red). 
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FIG. 5: Reweighting factor _R S from k s = 0.136 600 00 to 0.136 690 00 with Nb = 4 (black) and 8 (red) as a function of plaquette 
value. 
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FIG. 6: Reweighted plaquette value with i? s from n s — 0.136 600 00 to 0.136 690 00 as a function of the number of noise. Interval 
is divided by Nb = 4 (black) and 8 (red). 



1G 



K= 
P = 

- 'ss 
K- 



1 



II 



I 



10 



20 



t 



30 



0.8 



0.7 



0.6 



0.5 



0.4 



* 

A 

Z 
A 

N 



: - : 



I I I 



- -. : 



I I- 



10 



15 



t 



FIG. 7: Effective masses for the mesons (left) and the baryons (right) at the simulation point of (k u( j,k s 
(0.137 785 00,0.136 600 00). Horizontal bars represent the fit results with 1 standard deviation error band. 
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FIG. 8: Binsize dependence of the magnitude of error for the pseudoscalar meson masses. 
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FIG. 9: Effective masses for the mesons (left) and the baryons (right) at the target point of (fv* d ,K*) = 
(0.137 796 25,0.136 633 75). Horizontal bars represent the fit results with 1 standard deviation error band. 




FIG. 10: 7r, K and ?7 SS effective masses with the reweighting factors from (n u d,n s ) = (0.137 785 00,0.136 600 00) to (k^, kI) = 
(0.137 796 25,0.136 633 75). 
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FIG. 11: 7r, p and nucleon masses as a function of the number of noise. 
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FIG. 12: Hadron masses normalized by mo in comparison with experimental values. Target result for p meson locates below 
the figure. 
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FIG. 13: Determination of the physical point with ra,/mQ and mx/rofl inputs in (l/n u d, 1/k b ) plane. Solid and open black 
circles denote the original and target points, respectively. Green symbols represent (ft u d,tts) = (0.137 810 00,0.136 400 00) and 
(0.137 700 00,0.136 400 00) which are lightest two simulation points in Ref. [§]. 



